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COMPUTER PROGRAM FOR INVESTIGATING EFFECTS OF 
NONLINEAR SUSPENSION -SYSTEM ELASTIC PROPERTIES ON 
PARACHUTE INFLATION LOADS AND MOTIONS 

By Lamont R. Poole 
Langley Research Center 

SUMMARY 

A computer program has been developed to investigate the effects of nonlinear 
suspension-system elastic properties on -longitudinal vehicle and canopy motions and ten- 
sion loads experienced during the parachute inflation process. A mathematical elastic 
model of the suspension system is coupled to the planar equations of motion of a general 
vehicle and parachute canopy. Suspension-system elastic properties and histories of 
canopy drag parameter (product of drag coefficient and reference area), projected-area 
ratio (ratio of instantaneous projected area to projected area at full inflation), and added 
mass and its rate of change are input through the use of tables. The differential equa- 
tions of motion of the system are integrated by use of an equivalent fifth -order Runge- 
.Kutta technique. 


INTRODUCTION 

The successful operation of a parachute consists of two main phases: a deployment 
phase , in which the suspension lines and canopy are extended downstream of the towing 
vehicle, and an inflation phase, in which a net flow of air into the canopy alters the shape 
of the canopy until a final state of full inflation is reached. It is during this inflation 
phase that large forces are exerted on the towing vehicle by the parachute. 

It was shown in reference 1 that the elasticity of the parachute suspension lines can 
have an important effect on the loads exerted on the towing vehicle during the inflation 
process. In reference 1 the suspension lines were modeled as a linear massless spring 
with no viscous damping, and the analysis was limited to the case of infinite -mass (con- 
stant dynamic pressure) parachute deployment conditions. It was shown in reference 2 
that a massless -spring elasticity model provides a good approximation of the dynamic 
response of the suspension lines when loading conditions exist under which suspension- 
line wave mechanics can be ignored. Such conditions are normally in effect during para- 
chute inflation. The purpose of the present paper is to describe a computer program, 
PLAPIP (Planar Parachute Inflation Program), by which the analysis of parachute loads 



can be extended to include the effects of nonlinear spring characteristics, viscous damp- 
ing, and trajectory gradients on the loads transmitted to the towing vehicle. 

The function of program PL APIP is to solve the differential equations governing the 
planar motion of a general towing vehicle and the linear motion of an inflating parachute 
canopy relative to the vehicle. The two bodies are coupled by a two-component mathe- 
matical elastic model of the suspension system, which involves an additional equation 
governing the position of the juncture of the two suspension-system components. The 
differential equations are solved by use of an equivalent fifth -order Runge-Kutta integra- 
tion technique. 


SYMBOLS 

A reference area, meters^ 

C damping coefficient per member, newton-seconds 

Cj) drag coefficient, 

qA 

g acceleration due to gravity, meters/second^ 

h altitude, meters 

Kg ec specific secant modulus per member, newtons 

L unstressed length, meters 

m mass, kilograms 

n^ number of bridle legs 

n g j number of suspension lines 

q dynamic pressure, newtons/meter^ 

T s component, along the central axis, of sum of tensions in all suspension 

lines, newtons 

Tfo component, along the central axis, of sum of tensions in all bridle legs, 

newtons 
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t time, seconds 

v velocity, meters/second 

x distance measured relative to base of towing vehicle, positive rearward, 

meters 

/3 angle between central axis of parachute and any suspension line, degrees 

y vehicle flight -path angle, degrees 

£ . angle between central axis of parachute and any bridle leg, degrees 

p density, kilograms/meter^ 

Subscripts: 

a added (enclosed and apparent) 

bl bridle leg 

c canopy 

cs canopy skirt 

o initial conditions 

p confluence point 

si suspension line 

v vehicle 

00 free stream 

Dots over symbols denote time derivatives. 
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PROBLEM DESCRIPTION 


The problem considered in this paper is that of calculating the vehicle loads and 
motions which result from the inflation of a parachute attached to the vehicle by an elas- 
tic suspension system. It is assumed that histories of canopy projected-area ratio (ratio 
of projected area at any time to projected area at full inflation), drag parameter (product 
of drag coefficient and reference area), and added mass and its rate of change can be 
supplied from previous flight test data or auxiliary computer programs. The force 
transmitted to the towing vehicle is calculated by using elastic characteristics of the 
suspension-system components and a mathematical model of the geometry of the system. 
Vehicle trajectory and motion of the parachute canopy relative to the vehicle are 
determined. 

For the purpose of this analysis, the vehicle and canopy are considered to be mass 
particles, with the mass of the canopy particle varying because of changes in the air 
mass enclosed by the canopy during inflation and the "apparent” mass effect of changes 
in the kinetic energy of the air surrounding the canopy (ref. 3, pp. 153-155). The surface 
Of the planet is considered to be flat, and surface -relative accelerations are considered 
to be inertial. The motion of the towing vehicle is restricted to a vertical plane, and the 
motion of the canopy relative to the towing vehicle is restricted to occur along the central 
axis of the parachute. The central axis is considered to be parallel to the relative -wind 
velocity vector of the vehicle. 


Equations of Motion 

The forces affecting the motion of the towing vehicle and parachute canopy are 
shown in figure 1. 


The motion of the towing vehicle is influenced by aerodynamic drag, tension in the 
parachute bridle, and gravitational attraction. The acceleration of the vehicle in the 
direction of its velocity vector is given by 


v v = 


( c D A ) v q 


co + T b 


m. 


+ g sin y\ 


( 1 ) 


where q^ = | P^ 2 . 

In order to specify fully the surface -relative planar motion of the vehicle, the fol- 
lowing trajectory equations are required: 

h v = v v sin y 

y = ~g CQS-Z 
7 v v 
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( 2 ) 

(3) 



Vehicle velocity vector 



Figure 1.- Forces affecting motion of vehicle and canopy. 

The motion of the inflating canopy is influenced by aerodynamic drag force (CjjAq), 
tension in the parachute suspension lines, a force due to changes in the added mass, and 
gravitational attraction. The acceleration of the canopy in the direction of the vehicle 
velocity vector is given by 


■s - c q c - m a v c - m c g sin y 


m r + me 


where m a represents the added (enclosed and apparent) mass, and q c is assumed to 
be equal to | P^ 2 . 

In addition, the acceleration of the canopy relative to the vehicle can be written as 


x c = v v " v c 
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Mathematical Elastic Models of Suspension System 

In order to couple the equations of motion of the towing vehicle and inflating canopy, 
it is necessary to formulate a relationship between the tension in the bridle, the tension 
in the suspension lines, and the motion of the canopy relative to the vehicle. By con- 
sidering a particular suspension -system configuration and particular elastic characteris- 
tics, an equation can be formulated which defines the geometry of the suspension system 
and thus couples the equations of motion of the vehicle and canopy. 

For this analysis the suspension system will be considered as two sets of compo- 
nents: suspension lines and a multilegged bridle, which are joined at a common conflu- 
ence point. The suspension system is assumed to have negligible mass and negligible 
aerodynamic drag and is restricted to be under tension throughout the inflation process. 
The configuration for a fully inflated parachute, assuming axisymmetric inflation, is shown 
in figure 2. 



Figure 2.- Assumed configuration of suspension system. 

By using the geometry shown in the figure and assuming uniform elastic properties, 
an equation defining the component along the central axis of the sum of tensions in all the 
suspension lines can be written as a function of the suspension-line strain and strain rate 
as follows: 


T s = n sl cos^K^ 


( x cs " x p)sec 0 “ 

Lsl 


+ c sl 


(x c - Xp)sec p 


Lsl 


( 6 ) 


(Note that x c = x cs .) The secant modulus K g g ^ is a function of the strain, and the 
damping coefficient C s i can be a function of both the strain and the strain rate. 
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Analogously, an equation defining the component along the central axis of the sum of 
the tensions in the bridle legs can be written as 


T b = n w cos | 


/x^ sec | \ x 

KecJ\r- 1 ) +Cbl - 


Xp sec i; 


J bl 


(7) 


Under the assumption of negligible mass for the suspension system, the tension at 
any particular time is constant throughout the length of the suspension lines and bridle 
^T s = T^). Thus, by combining equations (6) and (7), an expression for calculating the 
rate of change of the relative position of the confluence point can be found, that is, 


*P 


— rc s i 5 

- l si l si 


Xp + K 

c sec 


,sl( ! 


cs 


Xp . Lsl cos ,)] - 

n sl C sl | n bl C bl 
L sl L bl 


( 8 ) 


This equation can be integrated to give the position of the confluence point relative to the 
towing vehicle. 

If suspension -system viscous damping is negligible (C s i = C b j = 0), equations (6) 
and (7) can be combined to form an algebraic equation defining the relative position of the 
confluence point, that is, 


X P = 


n sl K sec,sl(i^ - cos P) + n bl K Uc,bl cos * 

n sl K sec,sl r n bl K sec,bl 
L sl L bl 


(9) 


PROGRAM DESCRIPTION 


The entire computer program is written in FORTRAN IV for the Control Data 6000 
series computers. The main program, PLAPIP, is used for input and output of data and 
initialization of variables. Program PLAPIP uses Langley Research Center library sub- 
routine INT1A, which performs an equivalent fifth-order Runge-Kutta stepwise integra- 
tion of the system differential equations. (See appendix A.) Subroutine INT1A requires 
a user-supplied subroutine, CHSUB, which can be used for certain logical control. For 
the purposes of PLAPIP, CHSUB is a blank subroutine. Subroutine INT1A also uses sub- ' 
routine EQUOMO, which computes first derivatives of the system variables (v v , h v , y, 
x c , dx c /dt, and x p ^ at a particular time step based on tabular inputs of vehicle drag coef- . 
ficient, canopy projected-area ratio (ratio of the projected area at any time to the projected 


7 



area at full inflation), canopy drag parameter (product of drag coefficient and reference 
area), and added mass and its rate of change. Subroutine EQUOMO uses three library 
subroutines: AT62 furnishes data from the U.S. Standard Atmosphere, 1962 (appendix B), 
MTLUP performs first- and second-order interpolation from the tabular arrays 
(appendix C), and DISCOT performs a first-order (or higher, if necessary) interpolation 
to find intermediate values of a function of two independent variables (appendix D). 

Several options are available in the program. It will operate with either the U.S. 
Standard Atmosphere, 1962 (subroutine AT62) or user-supplied atmospheric tables. Also, 
the program allows for any values of viscous damping in the parachute suspension system. 

Main Program PLAPIP 

The main program, PLAPIP, is used for input and output of data, initialization of 
variables, and output of diagnostic messages. The flow diagram of program PLAPIP is 
as follows: 



Q 
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^ Print computed data 
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The program listing for main program PLAPIP is as follows: 


PROGRAM PLAPIP ( I NPUT , OUTPUT ♦ TAPE5=I NPUT * TAPE6=0UTPUT ) 

COMMON VAR (7) ,CUVAR(7> «DER (7) , N , NT , C I » Cl MAX , SPEC ,ELT , I TEXT ,ELE 1 (6) 
1 ,ELE2 (6) » ALTAB (50 ) ,RHOTAB(50) ,VSTAB(50) , NALT , XMTA6 (20 ) ,CDVTAB(20) * 
2NM,TINF(50) ,ARATIO(50> ,CDSTAB(S0) ,ADMASS(50) , ADMDOT < 50 ) » NT IN , STRSL 
3(50) , SLMOD ( 50 ) * NMOD * AMPSL ( 1 0 ) »SLRATE(10) ,NSLRT *CSLTAB ( 100) »NCSL,ST 
4RBL (50 ) »BLMOD ( 50 ) » NBLMOO » AMPBL (10) ♦BLRATE(IO) , NBLRT .CBLTAB ( 1 00 ) ,NC 
56L»WV* WCAN,G,OV*XLSL* XL6L»KEY'AT» TLONG, DMAX » XNSL * XNBL » ORAG, I SW* <3 IN, 
6QC,SLSTR»SLSTRR»XMIN*PA0ATP,XP,XPD0T»IPA, IPT.IPM* IPTL»IPT6 
DIMENSION TITLE (8) »ERRVAL (6) 

NAMEL I ST/T ABLE/ALT AB. RHOT AB, VSTA6,XMTA6,CD VTAB,T INF, ARAT 10 vCDSTAB, 
1 ADM ASS, ADMDOT , STRSL , SLMOD * AMPSL , SLR ATE, CSL TAB * STRBL « BLMOD * AMPBL « BL 
2RATE * CBLTAB ,NALT*NM, NT I N, NMOD » NSLRT » NCSL , NBLMOD , NBLRT ,NCBL 
3/1 NPUT/ WO, ALTO,GAMMAO*TIMEO»XCSO»XCDOTO»XPO»XPDOTO*G,DV»WV»WCAN»X 
4LSL,XLBL,DMAX,RADATP,XNSL,XNBL,ISW,KEYAT,TSTOP 
5/NINT/N,NT,CI,SPEC,CIMAX,ELEl»ELE2,ELT,ITEXT 
EXTERNAL EQUOMO,CHSUB 
PRINT 10 
10 FORMAT ( 1H1 ) 

20 READ(5, TABLE) 

IF (EOF, 5) 998,999 
998 STOP 
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999 CONTINUE 

WRITE (6 » T ABLE ) 

READ (5, INPUT) 

WRITE (6, INPUT) 

READ (5.NINT) 

WRITE ( 6 * NI NT ) 

READ 30 , T I TLE 
PRINT 40, TITLE 

30 EORMAT (8A10) 

40 EORMAT (1H1,8A10) 

PRINT 50 

C COLUMN TITLE FORMAT 

50 FORMAT (2X*TIME*4X*ALT*6X*VV*5X*VV-D0T*3X*M-INF*2X*Q-INF*3X*Q-CAN*4 
1X*XCS*5X*XCD0T*3X*STRSL*3X*STRRSL*5X*TL0NG*5X*DRAG*/2X*SEC*6X*M*7X 
2*MPS*6X*MPSS*10X*N/SM*4X*N/SM*6X*M*7X*MPS*5X*M/M*4X*PCNT/S*7X*N*8X 
3*N*/) 

C INITIALIZATION OF VARIABLES 
VAR ( 1 ) =T I MEO 
VAR (2) =VV0 
VAR ( 3) =ALT 0 
VAR(4)=GAMMA0*. 017453 
VAR ( 5 ) =XCS0 
VAR (6) =XCDOTO 
V AR ( 7 ) =XP0 

C INITIAL VALUES NEEDED FOR CALCULATION OF STRAIN AND STRAIN RATE IN EQUOMO 
XP=XP0 

XPDOT=XPDOTO 

C INITIAL VALUES OF PARAMETERS NEEDED BY MTLUP 
I P A=- 1 
I PT=- I 
IPM=-1 
I PTL=- I 
I PTB=- 1 
11=0 

60 CALL I NT I A (II ,N,NT,CI , SPEC » C I MAX » I ERR* VAR , CUVAR, DER » ELE l »ELE2 « ELT , 

1 ERRVAL .EQUOMO, CHSUB, I TEXT) 

C ANSWERS ARE ACCEPTABLE IF IERR= 1 OR 4 
GO TO (70,80,90,70) , IERR 

C CONVERT STRAIN RATE TO PERCENT-PER-SECOND 

70 STRATE=SLSTRR*1 00 . 

C PRINT-OUT FORMAT FOR COMPUTED DATA 

PRINT 75 , V AR ( 1 ) ,VAR(3> »VAR<2) ,DER(2) , XMIN , QI N , QC , VAR ( 5 ) »VAR(6> ,SLS 
1TR,STRATE,TL0NG,DRAG 

75 FORMAT ( 1XF5.3, 1 XF7 . 0 , 3XF5. 0 , 2XF8 .2 , 2XF5 . 2 ♦ 2XF6. 1.2XF6.1, 2XF6.2 * 3XF 
16.2,3XF5.4,2XF7.2,4XF6.0,3XF6.0) 

C CHECK TO SEE IF DESIRED TERMINATION TIME HAS BEEN REACHED 
IF(VAR(1) .LT.TSTOPIGO TO 60 
PRINT 76 

76 FORMAT (/2X*THE CASE HAS BEEN COMPLETED*) 

GO TO 20 

80 PRINT 81 

81 FORMAT </2X*ERR0R IN ELT BLOCK*) 

GO TO 20 

90 PRINT 91 

91 FORMAT (/2X*FAILURE TO ACHIEVE CONVERGENCE*) 

GO TO 20 

END 
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Subroutine EQUOMO 


Subroutine EQUOMO calculates values of the first derivatives of the system vari- 
ables at a particular time step. The flow diagram of subroutine EQUOMO is as follows: 


(equomo j 




* 
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Calculate 

array of 

derivatives : 

DER 

(2) = 

v v 

(eq. 

(D) 

DER 

(5) = 

hv 

(eq. 

(2)) 

DER 

(4) = 

Y 

(eq. 

(3)) .... 

DER 

(5). = 

* 

x c 


DER 

(&) = 

«* 

x c 

(eq. 

(5)) 


r-L- 

f RETURN 




The program listing for subroutine EQUOMO is as follows: 


SUBROUTINE EQUOMO 

COMMON VAR ( 7 ) *CUVAR(7) «DER<7> »N*NT*CI *CIMAX*SPEC*ELT» I TEXT »ELE1 (6) 
1 *ELE2 (6) ♦ALTAB(SO) »RHOTAB(50) *VSTAB<50) *NALT » XMTAB (20) *CDVTAB(20) ♦ 
2NM*TINF (50) »ARATIO(50) »CDSTAB(50) »ADMASS<50) » AOMDOT (50) *NTIN*STRSL 
3(50) »SLM0D(50) *NMOD*AMPSL (10) ,SLRATE(10) »NSLRT »CSLT AB ( 100 ) *NCSL*ST 
4RBL (50) *BLMOO(50) * NBLMOO ♦ AMPBL (10) ,BL RATE(IO) iNBLRT tCBLTAB ( 100)\NC 
5BL.WV*WCAN*G»DV*XLSL*XLBL.KEYAT»TL0NG,DMAX*XNSL*XNBL*DRAG»ISW*QIN* 
6QC.SLSTR*SLSTRR*XMIN*RADATP*XP,XPD0T» I PA* IPT* IPM* IPTL* 1PTB 
DIMENSION ANS ( 4) ,YT(4) *YA<2) ,VRT(.50*4) ,VRA(50*2) 

EQUIVALENCE (ARATIO* VRT) ♦ (RHOTAB* VRA) 

C*****KEYAT=0 DENOTES 1962 U S STANDARD ATMOSPHERE***** 

C*****KEYAT=1 DENOTES USER-SUPPLIED ATMOSPHERIC TABLES***** 

IF (KEYAT .EQ. 1 ) GO TO 10 

C THESE CONVERSIONS NEEDED BECAUSE AT62 USES BRITISH UNITS 
ALT=CUVAR(3)/.3048 
CALL AT62 ( ALT * ANS) 

RH0=ANS(l)*5l5.379 
VS=ANS(4)*.3048 
GO TO 15 

C INTERPOLATE IN USER-SUPPLIED ATMOSPHERIC TABLES 

10 CALL MTLUP(CUVAR(3) »YA«1*NALT *50*2*IPA«ALTAB*VRA) 

RHO=YA(l) 

VS=YA(2) 

15 QIN=.5*RH0*CUVAR(2)**2 

QC=.5*RH0* <CUVAR(2) -CUVAft(6) > **2 
XMIN=CUVAR(2) /VS 

C INTERPOLATE FOR AREA RATIO»CDS« ADDED MASS* AND M-DOT (YT<1) TO YT(4)) 
CALL MTLUP(CUVARd) *YT*2»NTIN,50*4,IPT.TINF*VRT) 
RAD=DMAX*SQRT(YT(1) )/3. 14159 
DRAG=QC*YT (2) 

CANMAS=WCAN/9.806 
TOTMAS=CANMAS*YT (3) 

C INTERPOLATE FOR VEHICLE DRAG COEFF. AS FUNCTION OF MACH NO. 

CALL MTLUP (XMIN*CDV»1»NM»20*1» IPM* XMTAB » CD VTAB) 

VDRAG=.7854*CDV*DV**2*QIN 

IF ( I SW . EQ. 1 ) GO TO 20 




C*****THIS BLOCK CALCULATES GEOMETRY » XP* AND TENSION FOR NO DAMPING**** 
XSL=CUVAR(5)-XP 
BETA=ATAN(RAD/XSL) 

REFLNG=XSL/C0S(8ETA) 

IF (REFLNG.LT. XLSUGO TO 25 
SLSTR=REFLNG/XLSL-1. 

SLSTPR=CUVAR(6) /XLSL/COS (BETA) 

XI=ATAN(RADATP/XP) 

BLSTR=XP/COS(XI J/XLBL-l. 

C INTERPOLATE FOR SPECIFIC SECANT MODULI 

CALL MTLUP(SLSTR,XKSL»2»NMOD,50»1»IPTL»STR5L»SLMOO) 

CALL MTLUP (BLSTR*XKBL«2»NBLM0D»50» 1. IPTB»STRBL»BLMOD> 

A 1 =XNSL*XKSL 
A2=XNBL*XKBL 
B1=A1/XLSL 
B2=A2/XLBL 

B3=A1*(CUVAR(5)/XLSL-C0S(BETA) ) 

C CALCULATE NEW POSITION OF CONFLUENCE POINT 
XP=(B3+A2*C0S(XI))/(B1+B2) 

TL0NG=A2*(XP/XLBL-C0S(XI) ) 

DER (7) =0. 

C*****END OF NO-DAMPING COMPUTATION BLOCK******************************* 

GO TO 30 

C*****THIS BLOCK CALCULATES GEOMETRY * XPDOT* AND TENSION FOR DAMPING**** 

20 XSL=CUVAR(5)-CUVAR(7> 

BETA=ATAN(RAD/XSL) 

REFLNG=XSL/COS (BETA) 

SLSTR=REFLNG/XLSL-I . 

IF (SLSTR.LT. 0. ) SLSTR=0. 

XI=ATAN(RADATP/CUVAR(7> ) 

BLSTR=CUVAR ( 7) /COS ( XI ) /XLBL-I . 

IF (BLSTR. LT. 0 . ) BLSTR=0 . 

SLSTRR= (CUVAR (6) -XPDOT) /XLSL/COS (BETA) 

BLSTRR=XPDOT/XLBL/COS(XI) 

C INTERPOLATE FOR SPECIFIC SECANT MODULI 

CALL MTLUP ( SLSTR ♦ XKSL » 2 » NMOD *50 » 1 . IPTL»STRSL» SLMOD) 

CALL MTLUP (BLSTR* XKBL » 2« NBLMOD ♦ 50 * I . IPTB * STRBL»BLMOD) 

C DOUBLE INTERPOLATION FOR DAMPING COEFFICIENTS VS. STRAIN AND STRAIN RATE 
CALL D I SCOT ( SLSTR .SLSTRR* AMPSL *CSLTAB , SLRATE* 1 1 •NCSLi NSLRT »CSL) 

CALL D I SCOT ( BLSTR .BLSTRR * AMPBL « C8LTAB » BL« ATE ♦ 1 1 »NCBL»NBLRT»C8L) 

A1=XNSL*CSL/XLSL 

A2=XNBL*CBL/XLBL 

A3=XNBL*XKBL*BLSTR*C0S (XI) 

A4=XNSL/XLSL 
A5=CSL*CUVAR (6) 

A6=XKSL* (CUVAR (5) -CUVAR (7) -XLSL*COS (BETA) ) 

C CALCULATE NEW VELOCITY OF CONFLUENCE POINT 
C**********0ER<7) IS EQUATION 8 ( XP-DOT ) ********** 
DER(7)=(A4*(A5*A6)-A3)/(A1*A2) 

XPDOT =DER ( 7 ) 

TLONG=COS (XI ) *XNBL* (XKBL*BLSTR*CBL*XPDOT/COS (XI ) /XLBL ) 

C*****END OF DAMPING COMPUTATION BLOCK********************************** 

I F ( TLONG. GT. 0 . ) GO TO 30 
25 TLONG=0 . 
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C CALCULATE DERIVATIVES 

C*********«OER (2) IS EQUATION 1 ( VV-DOT ) ********** 

30 DER(2)=-(VDRAG«-TLONG)*9.806/WV-G*SIN(CUVAR<4) ) 

C* ********* o £R ( 3 ) is EQUATION 2 (HV-DOT) *#**»**»»* 
DER(3)=CUVAR(2)*SIN(CUVAR<4> ) 

C**********DER(4) IS EQUATION 3 (GAMMA-DOT) ******* 
DER(4)=-G*C0S(CUVAR<4) ) /CUVAR(2) 

C*********»dER (5) IS RELATIVE VELOCITY OF CANOPY** 

DER(5)=CUVAR(6) 

VC=CUVAR (2) -CUVAR ( 6 ) 

C*******«**OER ( 6 ) IS EQUATION 5 ( XC-DOUBLE-DOT ) *** 

DER(6)=(DRAG+YT(4) *VC-TLONG*CANMAS*G*SIN (CUVAR (4) ) ) /TOTMAS+DER (2) 

RETURN 

END 


PROGRAM USAGE 


The program is run on the Control Data 6000 series computers under the SCOPE 3.0 
operating system and requires a field length of 40 OOOg storage locations. Most cases 
require a central processing unit (CPU) time of 50 seconds or less. However, an excep- 
tion to this rule is presented as sample case 1. This particular case involves a relatively 
long opening time and added mass effects and assumes no suspension-system damping. 


Input Description 

Input for the program is standard Control Data NAMELIST. Tables are listed under 
$TABLE; initial conditions, physical -system data, and switching parameters are listed 
under $INPUT; variables required by subroutine INTLA are listed under $NINT. Input 
and output are in the International System of Units. 

The $ TABLE input data are given as follows: 


ALTAB user-supplied array of altitude values h v for atmospheric data, m 


RHOTAB array of atmospheric density values p w corresponding to altitude values 
in ALTAB, kg/m^ 


VSTAB array of velocity of sound values corresponding to altitude values in 
ALTAB, m/sep . . , 


XMTAB array of values of Mach number (see CDVTAB) 


CDVTAB array of values of vehicle drag coefficient C^ v as a function of XMTAB 
TINF array of values of time, sec 
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ARATIO 

CDSTAB 

ADMASS 

ADMDOT 

STRSL 

SLMOD 

AMPSL 

SLRATE 

CSLTAB 

STRBL 

BLMOD 

AMPBL 

BLRATE 

CBLTAB 

NALT 

NM 


array of values of the ratio of instantaneous projected area of the canopy to 
projected area at full inflation as a function of TINF 

array of values of canopy drag parameter (Cj)A) c as a function of TINF, 

array of values of added (enclosed and apparent) mass m a as a function of 
TINF, kg 

array of values of time rate of change of added mass m a as a function of 
TINF, kg/sec 

array of values of strain in a suspension line, m/m 

array of values of specific secant modulus of a suspension line K gec 
as a function of STRSL, N 

array of values of amplitude of strain in a suspension line (see CSLTAB), 
m/m 

array of values of strain rate of a suspension line (see CSLTAB), (m/m)sec~* 

array of values of damping coefficient of a suspension line C s j as a func- 
tion of both AMPSL and SLRATE, N-sec 

array of values of strain in a bridle leg, m/m 

array of values of specific secant modulus of a bridle leg K gec ^ as a 
function of STRBL, N 

array of values of amplitude of strain in a bridle leg (see CBLTAB), m/m 

array of values of strain rate of a bridle leg (see CBLTAB), (m/m)sec~^ 

array of values of damping coefficient of a bridle leg C bl as a function of 
both AMPBL and BLRATE, N-sec 

number of values in ALTAB array (equal to number in RHOTAB and number 
in VSTAB) 

number of values in XMTAB (and CDVTAB) array 
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NTIN number of values in TINF array (also ARATIO, CDSTAB, ADMASS, 

ADMDOT arrays) - • 

NMOD number of values in STRSL (and SLMOD) array 

NSLRT number of values in SLRATE array 

NCSL number of values in CSLTAB array 

NBLMOD number of values in STRBL (and BLMOD) array 

NBLRT number of values in BLRATE array 

NCBL number of values in CBLTAB array * . - 

The $INPUT input data are given as follows: 

VVO initial velocity of towing vehicle, v v , m/sec 

ALTO initial altitude of towing vehicle, h v , m • 

GAMMAO initial flight -path angle of towing vehicle, y,-deg 
TIMEO initial time, sec 

XCSO initial displacement of canopy skirt relative to the towing vehicle, x cs , m 

XCDOTO initial velocity of canopy relative to towing vehicle, x c , m/sec 
XPO initial displacement of confluence point relative to towing vehicle, Xp 0 , m 

XPDOTO initial velocity of confluence point relative to towing vehicle,- x Pj0 > m/sec - 

G local acceleration due to gravity, g, m/sec^ 

DV maximum cross-sectional diameter of towing vehicle, m 

WV weight of towing vehicle, N 

WCAN weight of parachute canopy, N 
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XLSL 

XLBL 

DMAX 

RADATP 

XNSL 

XNBL 

ISW 


KEYAT 


TSTOP 

The 

N 

NT 

Cl 

SPEC 

CIMAX 

ELE1 


unstressed length of a suspension line, L s j, m 
unstressed length of a bridle leg, L^, m 
nominal diameter of parachute canopy, m 

radial distance from parachute central axis to point of attachment of bridle 
leg on towing vehicle, m 

number of suspension lines, n s j 

number of bridle legs, n^ 

parameter denoting viscous damping mode 

{ 0 no damping 

1 damping in suspension lines or bridle or both 

parameter denoting method of determining values of atmospheric conditions 

f 0 from U.S. Standard Atmosphere, 1962 

KEYAT = < 

[1 from user-supplied tables (ALTAB, RHOTAB, VSTAB) 
time at which user desires case to terminate, sec 
$NINT input data are given as follows: 

number of differential equations to be solved (for PLAPIP, N = 6) 
number of values in ELT array (for PLAPIP, NT = 1) 
initial computing interval for Runge-Kutta integration, sec 
time interval at which variable values are printed out, sec 
maximum computing interval allowed, sec 

upper bound of local relative truncation error for the respective 
dependent variables 
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ELE2 array of "relative zeros" for the respective dependent variables (absolute 

values below which relative error criteria are not applied) 

ELT a specific value of time at which control is to be returned to the main program 

(As PLAPIP does not require such a feature, an arbitrarily large number, 
which is greater than TSTOP, is input.) 

ITEXT integer code allowing for printout of history of computing interval and rea- 

sons for its variations 

ITEXT = 

In addition, a title card of up to 80 characters is inserted after the final card of the 
$NINT data. 

Output Description 

The output for program PLAPIP consists of first, the input data and then the com- 
puted data at time intervals specified by the input parameter SPEC. The input data are 
printed in the order in which they are read: first $TABLE and then $INPUT and $NINT. 
The title of the case is printed immediately preceding the column headings for the com- 
puted data. 

The computed data are presented columnwise as follows: 


TIME 

time, sec 

ALT 

altitude of vehicle, h v , m 

VV 

velocity of vehicle, v v , m/sec 

VV-DOT 

acceleration of vehicle, v v , m/sec^ 

M-INF 

free -stream Mach number 

Q-INF 

free -stream dynamic pressure, q^, N/m 

Q-CAN 

canopy dynamic pressure, q c , N/m^ 

XCS 

relative position of canopy skirt, x cs , m 



19 



XCDOT 

STRSL 

STRRSL 

TLONG 

DRAG 


Two sample cases are presented in order to illustrate the input and output quantities 
of program PLAPIP. The cases represent two inflations of a parachute, with 16.15-meter 
nominal diameter, in the Earth's atmosphere. One case is a subsonic inflation behind a 
slender towing vehicle at an altitude of 13.4 km, with the assumption of no suspension- 
system viscous damping. The second case is a supersonic inflation behind a bluff body at 
an altitude of 45 km, with suspension-system viscous damping included. 

Case 1 .- This case represents the inflation of a parachute, with 16.15-meter nomi- 
nal diameter, behind a towing vehicle, with 0.7 -meter diameter. There is no suspension- 
system viscous damping. Initial conditions are as follows: 

Mach number = 0.36 

h v = 13.4 km . 

q w = 1406.6 N/m2 

Y = -65.7° 

This case, having an opening time of 0.7 second, required a CPU time of 175 seconds. 


relative velocity of canopy, x c , m/sec 
suspension -line strain, m/m 
suspension-line strain rate, (m/m)sec~* 

component along central axis of total tension in suspension lines T s 
and bridle T^, N 

aerodynamic drag of parachute canopy, (CpAj c q c , N 

Sample Cases 


20 



The listing of the input data for case 1 is given as follows 

STABLE 
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THE CA5E HAS BEEN COMPLETED 



Case 2 .- This case represents the inflation of the same parachute behind a towing 
vehicle with a 3.505-meter diameter. The suspension-line viscous damping coefficient 
is a constant equal to 100 N-sec per line. Initial conditions are as follows: 

Mach number =2.16 
h v = 45 km 
q M = 488.6 N/m 2 
y = 16.0° 

This case, having an opening time of 0.185 second, required a CPU time of 17 seconds. 
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The listing of the input data for case 2 is given as follows: 
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The sample output for case 2 is given as follows 


w 

— X 

■x 

<—> 

10 

A 

r-H 

OJ 

4 

X 

4 


n 

IO 

in 

A 

X 

OJ 

1 ^ 

to 

X 

OJ 

4 

X 

to 

4- 

X 

<y 


O' 

X 

iO 

in 

r-f 

4. 

.in 

r> 

4 

X 


X 


X 

1 -^ 

C^ 

X 

X 

X 

r- 

ro 

IO 

4 

X 

O 

10 

A 

O 

X 

ro 

X 

4 

X 

O 

0 


«-H 

O' 

X 

a ^ 


»—* 

A 

ro 

iO 

X 

X 

O 

ro 

IO 

X 

OJ 

IO 

0 


Si 

IO 

X 

(VI 

X 

04 

X 

X 

ro 

O 

X 

X 

to 


O' 

X 

to 

4 

0 








<—* 

»— 

r— 

—+ 


r— < 

0 J 

A 

ro 

4 

IO 

h- 

X 

O 


ro 

4 

X 

X 

OJ 

0 - 

ro 

O' 


r- 

X 























i— 1 




1 -^ 

OJ 

ru 

ro 

ro 

4 

to 

X 


X 

0 O — h- r- — • A 

G 

r^- 

X 


*— < 

0 

4 

0 

r-t 

O' 

O' 

G 

O' 

r~ 

r-^ 

r- 


0 

OJ 

X 

X 


r- 

S 

X 

0 

(VJ 

s 

ro in 

X 


X 

OJ 

O' 


4 

ro 

4 

X 

X 

4 


A 

4 

1 ^ 

4 

0 

«— 

X 

o- 

r- 

X 

ro 

o- 

r-t 

ru 

0 e 



r^H 

r— 

(V 

OJ 

ro 

4 

to 

X 

r- 

O' 

0J 

10 

O 

X 

10 

X 

O' 

r— 

X 

X 

X 

4 

<—* 

X 

O' 

r- 

-1 















rj 

A 

ro 

4 

*x 

O 

ro 

r- 

i-~* 

4 

rc 

(\i 

X 


K 




















f— 

*— 1 

— < 

(VI 

OJ 

OJ 

ro 

ro 

4 


_J x 

in x 
a h- 
cc z: 
o 

(/> X 


o x a O' x a o 4 o c a x ro a x g x 4 ^ 0'4J' O'X — ^ r- xg£ x a 
CpH(\j^jMrt^4 , a^(Vj-HOMr'C( r )^i?OC^fV;C{V!OCnvC(VC^nfnvtC 


I I I I 


^(\jf r )4 ->tJi«£)Q 0 onNMccin^Nninx 4 a ro p- x 
-< h h !\j i\i n 4- x x to x> x >- -c o' o a 


_j 0 r— r-<A)X 4 X£r'- 0 '-— roxo 4 acxo' o 4 x g o t' a O' •— x x ^oo r- y 

(/) r oooogogooo-« — I H ,M (\J N n P) ^ T -C T T H .} T O n N H nC (\J t .T 

Ll \ oggooogoooooooogoooogog^— < -• (\, ..\j f\j .'n n 4 4 xi 

►— 7" ooooooooooooooooooooooooooogoogooo 

x • • • • 


G 

O h~ 

— O X 
c a 

II O .7 

X 


o ^ £^n ( \j^4OHf)N{\iHiNooMa v OL0aH^!\)'C'^irf)T'C^ p- 
gogggo — aj 4 x x — 4 x x a x : x 4 r, c vC 4 -4 -« r< 4 it rn "0 -4 


7 — I I I I 


— . — . — — A A 


4 TN t crt!4£NTO-r.^Cif 
— « — • — • tv- (\i iv n r, 


Jj 

u 

>— X 
X u 


x x t> x x x x x x- v O' O' o — a co 4 it r- o a r- Ai p 4 \ o y y .y g ai x x 

o c o o g c. - c. g g g g — « —• »-i -« — ^ A a- svi x X 4 X« X X p- x o — ■ rvj 4 

C'OGGriOCC-GGGGGOGOGO'rGOOGOGGOGGG-— — < — — - 

a. r rp rn a ro ^ r p x a x> to p p, p p p p p p p r p p a- x< r ^ P n a a 


S < 7 
»- U X 
A I X 

j 3 ^ 
< 


X X 4 A — 7 £ X T' X 3 4 O' o PT PT fV - > C ^ P X.T M 4 ■H N 4 X- -X 

X X 4 ' a. a. N P > £■ £ X X 4 4 4 A (V — O X XX,- c X (V c: J N -X — A — 4 

X rc x x x. x x x x x x x x x x x x x r- r- r- p- p- x x x: x m x, x. r 4 <4 a 

4 44-4 4444 44 444 -4 444 44 4 44 4 4444444444 4 


-c ;- 

C — IT 
~ I X 

a o Z 


•x x 4 a — 7 x* x x a (vs o O' x r x. n m o ■> n x, 4 a g x x a o p- a o x o 

•X X X -X t P N P P p* P N vC ^ i' 'f X XX XX X X ‘.O X 4 4 4 -4 ‘X A (Vi A A 

X X X ao X X X X X X X X X T X X X XX T> X X X X r X X: X X> X X £ X X 

44 4 5 4 4 4 4444 44 4 444 4 4 4 44 4-4 4444 4 4 4 4 44 


x x £■ x x x x x- £ x x x x x x- r x c c • x c x x x x c c x r x x c r x 


■J 

1 

• 

VI 

• 

:\J 

• 

V 

• 

■\J 

« 

rvj 

• 

At 

• 

• 

■Vi 

• 

rvi 

• 

'\J 

• 

AJ 

• 

OJ 

• 

AJ 

• 

AJ 

* 

AJ 

* 

A 

• 

V' 

• • 

■\; a 

• 

■V) 

• 

\J 

• ■ 

V A' 

• 

AJ 

• 

AJ 

• 

AJ 

• 

A 

• 

AJ 

• 

A 

• • 

A AJ 

• 

A 

• 

A 

• 

A! 

X 

1 - X 

o- 


X 

r- 

X 

y 

. 

■t 

X 


T 

r^ 

X 

X 

J 

r» 

J 

P- X 

■/ 

r- 

X' o* 

t 

Y*‘ 

3 ' 

A 


S 


X 

ro 

V, 

r 

7 X 

5 

.-r 

1 

5 

i 

t 

r 

r 

X 

c 

X 

P> 

X 

■> 

G 

— . 

O 

T 1 

•A; 


1 0 

X 

A; 

— 

A! 


>.r 

7 V 


J 


u.’ 

r_‘ ^ 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

■ 

• • 

• 

« 

• t 

• 

• 

• 

< 

• 

• 

• • 

■ 

« 

• 

a 

1 7 

. — ■ 

. — 


r— » 

— • 

— < 

— » 

— . 

— 1 

— ■ 

p— 



— * 

A 

'A 1 

;V 

0 : A 

P" 

ro 

1 r 

X 

•X 



•A' 

X 

7 7 

4 

7 

w 

■ 7 . 

> 

— 

— 

— 

— * 

•— 

r-H 

— • 

— 

— * 

— 

■— 

—4 


— 

— < 

— 

— 

<— ■— 


— 

— • — 

— • 

— • 

■A: 

A' 



*0 1 

4 

X 

X 


cr ■> 

i 


i i i i i i 


i i i i i i i i i i i i i t i i 


i i » i 


h~ ir- 
'S > y. 


a 


■r r. r r it s ip t / 444444 1 4 4 4 4 4 4 4 4 r a ^ a ^ ^ ~ a p v 

r — r — , — r — p — p — fs. 1^ r^. ^ r— r — r — r— r^-. r— p~ 





cr 

r— 

.A: 

r~- 

7 

P 

X 

r^. 

A 

■7 

O 

— i 

A: 

p 

4 

•0 

:. 

X 


■1 

T 

'P 

— 1 

A 

- y *- 

4 

O 

£ 


X 

•t n 

— 

A 

■O' 



0 

C' 

O 

cr. 

O 

c~ 

O 

O 

cr 

O 

r— . 

— . 

— 1 

^-1 

— < 


—4 

r-H 

«— 1 

— < 

— < 

A 

a: 

A 

A 

P 

A- 

A 

•A, 

;\ 

>\ A 

A) 

A 

u 

; 

£• 

c 

cr 

C 

G 

O 

c 

_■ 


0 

O 

C 

0 

rz 

c; 

s> 

c 

O 

£; 

cr 

• “ 

c_ 

G 

0 

G 

- 

G 

G 

£ 

CO 

£. 

G G- 

.“• 



<*■ 


X 

!X 

P 

X 

X 

•0 

X 

X 

X 

10 

X 

X 

X 

X 

to 

0 

X 

X 

X 

X 

0 

O 

X 

X 

X 

O 

X 

O 

X 

X 

O X. 

X. 

O 

1 



4 

4 


t 

4 

7 

4 

t 

■t 

4 


4 

4 

4 

t 

7 

4 

4 

4 

4 

4 

..*» 

4 

t 


4 

4 

4 

4 

'7 

5 4 

4 

4 

: f 

UL-' 


O 

X 

O 

X 

O 

X 

O 

tr. 


X 


X 

0 

X 

0 

X 

_. 

X: 

O 

X 

O 

lO 

O 

ir 


X 

O 

X 

0 

X 

G X: 

0 

X 

• 

2 

(J 

0 

O 



•A 

A 

r~ 

ro 

4 

4 

to 

X 

X 

X 

p- 

p- 

or 

a 

7 

O' 

G 

C 

. — ■ 


A 

A 

P' 

A 

4 

4 

X X 

£ 

£ 

X 

1 — 

LlJ 

O 

O 

O 

0 

O 

O 

0 

0 

O 

0 

0 

O 

O 

0 

0 

0 

0 

0 

O 

O 

— * 

r— 

*— 

. — 

r-* 


r-f 

— 1 

— < 


r— t — < 

r— 

r— 1 

■— 1 

t— 

X 

« 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

* 

■ 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• _ • 

• 

• 


30 



170 4S0.33 701 -66.39 d.lo 431. 6 4?6.1 Jl.6,2 4 1 . 63 .0041 14H.S0 4631 7 76252 


O O' O' C C ? 00 ip O' J MT, n M 

rv a m cv o' r- <x x -h p- x r- 4- x o .x> n 

<r®(\jNnO(\JNNv£)ajXCC 4 n(VJ 
m x a ^ohhc^^d^ctnncvC 


co h x, ,0 tcnNnM>oxNO'oji 
LT S a f\' C 4 ■‘Too^xooxa'oxm 

»-• 4 a? it. -h r- oxxo4xx— <xa'4 
^ x o n co f\j \0 o ' — 1 o rvj 4 in 
4XX'Cf^f^r^aoaoa'0'a' x cc p- x x 


cr rv x x 4 r-* rv'oo'oo^vtar^ox 
1/1 !V n n o n ^ t m 4 - iX '-*— ^ 1 r- uo — * 

4 o x m cr HrtMXNH^njioir 

n H ro M 'J 1 X iT, H £ c 4 :\f O' X — D 
-^r\jf\jf\if\jrvif\jnjoj^H^ 1 1 . — « rvj r\j 

1 1 1 


4 > x — * <m -x x r- x x x <x rvo c 4 no 

m 4 X 0 '-* 4 ^O'X'X- 4 r- X X X O' r^ 
r-xo' 0 (vm 4 LDX( s - xxxx-^xd 


o' d — •^ffco^mnifoif, x o r\j 4 

r-x 4 rvj O' d r- n 4 O' x <vj rvj x 

■x -x a x o — • x r vj — t x O' — • x r- 4 o nj 

4 iD O i'*- X X P* r- x 4 PJ — < 1 4 X r- 

1 1 1 1 


LD ro 0 ' 0 ' o O 00 — ‘ O' X X O' C OJ x aj 
x ^ 4 N* ^ X or r )f v >j'<— <(\j(\i{\ior s ~ 4 


• w <\j M n m 4 44 4 in it m :D ir 4 4 

XX P p P ^ p p .r i»; p r, ^ (X 


4 4 MT M 7 1 XX 4 0* O i 


4 0 4 X O' 


x m o o -4 (Vi D rv it cvj 4 x mo'm 4 — 
— oo'Xf^r^r^xO'^roxJcooroiDN 
4 4 m m m mm m <m 4 44 4 x x x x 


o tx 0- rvj in x 


m id x O' — • <m id a* oj 


— I O a 0 * X N N X' l/i P Ai M ^ c o o 
444444444444444 44 


•D D X D D 4 4 4 4 4 4 m m m m . 


oj xj rvi rvi xj x> rvi <\j xj ,v xj xj ,xj xj xj xj v» 


t o -i n o m o d x r- X — • r- x torn 

c x xi n — « m 4 x — -xj y o 40 4 — • — • 

,j; 4 — • d m p*» cv/ r- rv x r* p- ^ p*» x x. x ; 

x^x r o o o c — * —1 ^ — « o o a h» 

| | | | | l | | | 

I I t I I I I I 


— — • O e O O- T X h- r- X X X T 4 4 c: 

c-ooooooo'oocrj'oa'ooj 
x x x x x x x x x x x x x 


4 if! x r- x ^ o o — • cv m 4 x x n- a. a- 

r. m r*. m m 4 4 - 4 - 4 4 * -4 4 -4 4 4 4 

00000000 00 o o c o c © o ■ 

ID ID X X ID X ID X X X X, X 'D X X X X 

444444 4 4 4 4 44444 4 4 


X OX OX OX OX OXO XOX OX 

r- cc a a cr 0 o ^ ^ (Vi fv m m 4 4 x x 
^^-.^^H^Ojr^ojnjojnjrvjDinjoj (\i(\j 


♦pH 

I • 

< £ 
S2 

a ~ 

rv O 

ra M 


u,' 

1— 

Ui 

_J 


d 

d 

X 

c 

d 

a 

2 : 

a; 

0 

■ *-H 

-M 

d 

~ 


d 

> 

u 

d 

oi 



<D 

U 

d 

0 

d" 

0 

Id 


*-< 

-4-> 

X 

x: 

0) 

a 

X 

0 

Sh 

C 

a 

l r 

d 

pH 

rt 

< 

<D 

a* 

K 

X 

Cfl 

d 


CD 

0 


X' 

X 

PS 

•pH 

~8 


< 

0 

>> 

<u 



pH 



X 

b£> 



X 

d 




31 



APPENDIX A 


LANGLEY LIBRARY SUBROUTINE INTI A 

Language : FORTRAN 

Purpose : ENT1A is a closed subroutine for the solution of a set of ordinary differential equations. 

Use : CALL INTI A (II, N, NT, Cl, SPEC, CIMAX, IERR, VAR, CUVAR, DER, ELE1, ELE2, ELT, 

ERRVAL, DERSUB, CHSUB, ITEXT) 

II INTI A is composed of an initialization section and an integration section. The user is required 

to enter the initialization section before he starts his first integration step. The above calling 
sequence is used for both initialization and integration with the value of the code word n deter- 
mining which of the two sections of INT1A will be entered. 

The user must set II = 0 in order to initialize. 

During initialization the derivatives will be evaluated using the initial values of the variables but 
no integration will occur and control will be returned to the calling program. When INTI A is 
called with n > 0, entry is made to the integration section. Upon each entry to INTI A, the sub- 
routine stores a 1 in II so that the users need not supply a value of II > 0 for repetitive 
integration. 

Besides serving as a means for specifying the entry point to INTI A from the calling program, II 
can also be set to specified values in CHSUB to accomplish the following: 

2 The user will store the integer 2 in n if the answers in CHSUB are not acceptable to him 

and he wishes to recompute the answers using a shorter interval. This shorter inter- 
val must be stored by the user in CI. It must be smaller than the computing interval 
just used. 

3 The user will store the integer 3 in n if he wishes to return to the calling program. 

The answers for the interval are considered acceptable to the user and will be trans- 
ferred to the VAR array (explained below) by INTI A. 

In DERSUB, II may be set to: 

4 The user will store the integer 4 in n if he wishes to discontinue calculation of the pres- 

ent interval and return to the calling program. On return to the calling program, the 
answers at the beginning of the interval will still be in the VAR array. 

If the user does not set H to a value in either CHSUB or DERSUB, H will always be 1 upon the 
return to the calling program. 

N An integer value supplied by the user which is the number of differential equations to be solved. 

Subroutine INTI A is compiled to solve a maximum of 20 equations but may be recompiled for 
larger values of N if necessary. 

NT An integer value supplied by the user which is the number of values in the ELT block 

described below. Subroutine INT1A is compiled with a maximum of 10 values in the ELT 
block but may be recompiled for more values if necessary. 
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APPENDIX A — Continued 

Cl A floating-point value supplied by the user which is the computing interval INTI A will use ini- 
tially. Cl must be a signed value, positive if integrating forward, negative if integrating back- 
wards. Upon entry to CHSUB, Cl will contain the computing interval that INT1A will use for 
the next step unless it has to take a short interval to hit an ELT value or a SPEC value 
described below. The computing interval used on the present step is available in CHSUB as 
the algebraic difference between CUVAR(l) and VAR(l). Since the subroutine is used on a 
binary computer and the interval variation is a halving and doubling process, Cl should be a 
power of 2. 

SPEC A floating-point value supplied by the user which specifies how often he wishes INTI A to return 
control to the calling program so that the user may print his results. 

SPEC = 0.0 Control will be returned after every acceptable integration step. 

SPEC > 0.0 SPEC is the absolute value of the specified increment of the independent variable for 
which the user desires control returned to the calling program. 

The first printout is made at the initial value of the independent variable. The next return is at 
the nonzero integer multiple of SPEC closest to the initial value of the independent variable. The 
remaining returns occur at values which have been updated from this point by the increment given 
in SPEC. The return times generated by the increment given in SPEC are not altered by an inter- 
vening return due to an ELT value (explained below). 

CIMAX A floating-point value supplied by the user which is the absolute value of the maximum computing 
interval that will be used. This value will be used if the doubling process would extend the 
computing interval to a value larger than CIMAX. CIMAX should be set to 0.0 if there is no 
desired maximum. 

IERR An integer value supplied by INTI A as an error code. It must be checked at every return to the 
calling program. It may have the following values: 

1 A normal return, no error. 

2 The ELT block is not monotonic in the direction of integration. 

3 The variables have failed to meet the local truncation error requirements nine 

consecutive times. The answers at the beginning of the interval are still in the 
VAR array. 

4. The variables have failed to meet the local truncation error requirements at 

least nine times over the last three intervals. An acceptable answer has been 
reached, however, and is in the VAR array. 

VAR A one -dimensional array containing the independent variable followed by the N dependent vari- 
ables. The user must store the N+l initial values in the array for initialization. INT1A will 
store the new values of the variables in VAR after each integration step when they are 
accepted by the user in CHSUB. The elements of the VAR block can be printed out in the call- 
ing program in accordance with the user's specification in SPEC. 
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APPENDIX A - Continued 


CUVAR 


DER 


ELE1 


ELE2 


ELT 


ERRVAL 


DERSUB 


A one -dimensional array which is given values by INTI A for two purposes. INTI A will store in 
the same order as the VAR array the values of the independent variable and N dependent vari- 
ables at which it wishes the derivatives to be evaluated in the DERSUB subroutine. 

INT1A will also store the tentative answers after each integration in the CUVAR array before 
calling CHSUB so that the user can check these values to decide to accept or reject the 
answers. If accepted, the CUVAR values will then be transferred to the VAR array. 

No values need to be initially stored in CUVAR. 

An N+l array in which the user will store the derivatives evaluated in DERSUB. The derivatives 
should be arranged by the user in DERSUB in the same order as the VAR block so that DER(2) 
will be the derivative of the variable stored in VAR(2), and so forth. DER(l) will be unused. 
The derivatives must be computed using values of the variables which have been stored in 
CUVAR (not VAR) in INTI A. 

A one -dimensional array of N values supplied by the user each of which is the upper bound of 
local relative truncation error for the respective dependent variables. If the error for any 
variable exceeds its respective ELE1 value, the computing interval is halved and the integra- 
tion restarted at the beginning of the present interval. If the error for all of the variables is 
less than 1/128 of their respective ELE1 values, the computing interval is doubled for the next 
integration step. 

A one -dimensional array of N values supplied by the user which represents a small value of 
"relative zero" for the respective dependent variables. If the absolute value of any of the 
variables is less than its respective ELE2 value, the relative error criteria for that variable 
will not be applied. 

A one -dimensional array of NT values supplied by the user which are values of the independent 
variable at which the user specifically desires control returned to his program. The values in 
the ELT block must be monotonic in the direction of integration or an error return will be 
given by INTI A. 

A one -dimensional array of N elements in which INTI A stores an estimate of the local truncation 
error for each of the N dependent variables. The relative errors are computed from these 
values and compared with the specified ELE1 values. 

The name of a subroutine written by the user which will be called by INTI A to evaluate the 
derivatives. The derivatives must be stored in the DER array. INTI A will call DERSUB to 
evaluate the derivatives with the values of the variable it has stored in the CUVAR array. 

The name given to the DERSUB subroutine must appear in an EXTERNAL statement in the call- 
ing program. The user may return to the calling program by setting II to 4. 
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APPENDIX A - Concluded 


CHSUB The name of a subroutine written by the user to allow certain logical control. After each integra- 
tion step, INTI A will make available to the user in CHSUB the tentative answers in the CUVAR 
array. The VAR array will contain the last accepted answer (that is, the value of the variables 
at the beginning of the interval). Whenever the user specifies the answers are acceptable, the 
values in the CUVAR block are transferred to the VAR block. In CHSUB the DER block will 
contain the values of the derivatives evaluated with the present CUVAR block. The user has 
three options: 

(1) Not change II. II = 1 is considered by INT1A to denote that the user has accepted the 

answers in the CUVAR block. II always equals 1 upon entry to CHSUB from INT1A. 

(2) Set II = 2. The user does not accept the answers and wishes to recompute the interval 

using a new computing interval which he stores in Cl. This computing interval must 
be smaller than the computing interval just used. This new value of Cl will now be 
stored by INTI A as the normal computing interval for the subsequent integration steps. 

(3) Set II = 3. The user accepts the answer but wishes to denote a condition that he can test 

in the calling program. Control will be returned to the calling program with the 
answers in the CUVAR array transferred to the VAR array. 

The name given to the CHSUB subroutine must appear in an EXTERNAL statement in the call- 
ing program. 

ITEXT An integer code word supplied by the user which gives him the option to have INTI A print out a 
time history of the computing interval and the reasons for its variation. This printout should 
be requested only for problems which must be rerun because of unsatisfactory results the 
first time. 

ITEXT = 

Restrictions : See arguments listed under CALL statement. 

Method: Subroutine INT1A, written in coordination with the other integration subroutines in the INT(x) 
common -usage series, is a fifth -order integration subroutine. The classical fourth -order Runge- 
Kutta formula is applied in conjunction with Richardson's extrapolation to the limit theory. INT1A 
is a variable interval size routine, in which the interval is varied to meet a specified local relative 
truncation error. 

Accuracy : The variable interval size mode of logic is used to make available an estimate of the local rela 
tive truncation error which is then controlled as explained in the ELE1 block discussion. 

Storage : 2557g locations. 

Subroutine date : August 1, 1968. 


1 0 No printout requested 
\l A printout requested 
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APPENDIX B 


LANGLEY LIBRARY SUBROUTINE AT62 

Language : FORTRAN 

Purpose : Subroutine AT62 approximates the U.S. Standard Atmosphere, 1962. It computes density in 
slugs/ft 3 , pressure in lb/ft 2 , temperature in degrees Kelvin, and the velocity of sound in ft /sec at any 
geometric altitude in the range between -16 500 feet and 2 320 000 feet 0 

Use : CALL AT62 (Z,ANS) 

Z Geometric altitude in feet 

ANS A one -dimensional array that contains the results: 

ANS(l) Density in slugs/ft 3 

ANS(2) Pressure in lb/ft 2 

ANS (3) Temperature in degrees Kelvin 

ANS (4) Velocity of sound in ft/sec 

Restrictions : For altitudes below -16 500 feet the values of density, pressure, temperature, and velocity of 
sound are not valid. The concept of the velocity of sound in the atmosphere becomes essentially meaning- 
less at altitudes in excess of 300 000 feet. To point out this limitation, the velocity of sound at altitudes 
above 300 000 feet is set equal to the velocity of sound at 300 000 feet. For altitudes above 2 320 000 feet, 
density, pressure, and temperature are set equal to their respective values at 2 320 000 feet. 

Method : The equations and techniques are identical to those used in computing the U.S. Standard 
Atmosphere, 1962. 

Accuracy : The tables in reference (a) were computed with an IBM 7094 by use of some double -precision 
arithmetic. In converting the routine for the Control Data 6000 series computers, all double -precision 
arithmetic was eliminated. Accordingly, there may be slight differences between the results of the con- 
verted subroutines and the tables. 


Refer ences : (a) Anon.: U.S. Standard Atmosphere, 1962. NASA, U.S. Air Force, and U.S. Weather Bur., 
_ ” Dec. 1962. 

Storage : 16540 locations. 

Subroutine date: August 1, 1968. 
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APPENDIX C 


LANGLEY LIBRARY SUBROUTINE MTLUP 

Language : FORTRAN 

Purpose : Multiple Table Lookup (MTLUP) computes Yj = Fj(X) for j = 1, 2, . . NTAB from a set of 
NTAB tables using first- or second-order interpolation. An option to give Yj a constant value for any X 
is also provided. 

Use : CALL MTLUP (X, Y, M, N, MAX, NTAB, IP, VARI, VARD) 

X The name of the independent variable X 

Y The name of a one -dimensional array of the dependent variables calculated, Yj = Fj(X) for 

j = l,2, . . NTAB. 

M Interpolation parameter (integer), 1, 2 for first- or second-order interpolation, 0 for Yj a 

constant as explained in the section entitled "Notes." 

N The number of points in each of the tables (integer), N ^ MAX. 

MAX The maximum number of points in each of the tables (integer) as given in the DIMENSION state- 

ment in the calling routine. 

NTAB The number of tables of dependent variables (integer). 

IP Interval pointer (integer variable name, not a literal integer constant), unique for a given VARI, 

with the following three uses: 

(1) Prior to interpolation. - If IP = 1 prior to interpolation, the values in the VARI array are 

checked to determine if the entire table is strictly increasing or strictly decreasing. 

Note that equal values in the VARI array will give an error condition. 

(2) For interpolation.- Upon entry, 1 = IP £ N - 1 gives the interval (VARI(IP), VARI (IP + 1)) 

to be checked first for X. 

(3) After interpolation.- Upon return, 0 i IP 5 N gives the interval where X was found; 

VARI(IP) £ X £ VARI(IP + 1) for strictly increasing, VARI(IP) £ X ^ VARI(IP + 1) for 
strictly decreasing. IP = 0 indicates low-end extrapolation and IP = N indicates high- 
end extrapolation. 

If MTLUP is used for more than one independent variable table (VARI) in a single program, a 
different integer variable name should be given for IP for each VARI array used. Each 
should be initialized to -1 and then left in variable form. (See section entitled "Method.") 

VARI The name of a one -dimensional array which contains the N values of the independent -variable (X) 
table. 

VARD The name of a two-dimensional array in which each of the NTAB columns contains the N values 
of a different dependent -variable (Yj) table. 
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Notes: VARI(I) corresponds to VARD(I,1), VARD(I,2), . . and VARD(I, NTAB) for 1=1,2,. . N. For 
M = 0 or N i 1, Y(J) = VARD(1,J) for J = 1, 2, . . ., NTAB, for any value of X. The program extrapolates. 

Restrictions: The following arrays must be dimensioned by the calling program as indicated: VARI(MAX), 
VARD(MAX,NTAB), Y(NTAB). 

The independent variable X and all the values in VARI and VARD must be floating point. M, N, MAX, 
NTAB, and IP must be integers. The values of the independent variable X in the VARI array must be 
strictly increasing or decreasing. 

Method: The use of N i MAX allows the user to specify a maximum dimension (MAX) for the arrays and 
then run cases with smaller actual dimension (N) without recompiling, by reading N as program input. 

The interval. pointer IP is used to check the monotonic order of a given table only once. Thereafter, it 
preserves the location in the table of the previous answer. 

In table lookup for scientific computing, consecutive X arguments tend to fall in the same region of 
independent -variable table. To take advantage of this, IP provides, for each VARI, a storage location 
to be used by the subroutine in communicating from one call to the next. When X is found, IP is set to 
the interval count where X is located. Then, on a subsequent call with the same variable name for EP 
(and, consequently, the same VARI table), the interval search begins with the interval that contained the 
previous X. This feature generally saves a significant amount of run time. 

If x i,X 2 , . . ., X N are the tabulated values of the independent variable and Y l5 Y 2 , . . ., Y N are the 
corresponding values of the dependent variable, then the interpolation equations (derived from refer- 
ences (a) and (b)) are as follows: 

Second order: For i chosen such that | X - Xf | + |X - X^ +2 | is a minimum, 
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APPENDIX C - Concluded 


First order: For Xj i X § X 1+ j , 

(T M -Y,)(X-X|) 

' Xitl - X, 

Example: Given a table of 30 values of altitude (ALT) and tables of 30 values each of 
temperature (TEMP) and velocity (VEL) as functions of altitude, 

DIMENSION TEMP(50), VEL(50), Y(2), ALT(50), VARD(50,2) 
EQUIVALENCE (TEMP, VARD(1,1)), (VEL, VARD(1,2)) 

IP = -1 
X= 302.6 


CALL MTLUP (X, Y, 1, 30, 50, 2, IP, ALT, VARD) 
PRINT 10, Y 


X = 198.7 

CALL MTLUP (X, Y, 2, 30, 50, 2, IP, ALT, VARD) 
PRINT 10, Y 


Accuracy : Accuracy is a function of the order of interpolation (M) used. 

References : (a) Nielsen, Kaj L.: Methods in Numerical Analysis. Macmillan Co., c.1964, pp. 87-91. 

(b) Milne, William Edmund: Numerical Calculus. Princeton Univ. Press, 1949, pp. 69-73. 

Storage : 444g locations. 

Error condition : The first call to the subroutine using a new VARI table (with IP = -1) checks the monotonic 
order of that table before interpolating. If an error in the order is found, the subroutine will print 
TABLE BELOW OUT OF ORDER FOR MTLUP AT POSITION xxx X TABLE IS STORED IN LOCATION 
xxxxxx (absolute). It then prints the contents of VARI, and STOPs the program. A call with 0 ^ IP £ N 
will interpolate but will not check the monotonic order of VARI. 

If N = M = 2, the M TABLE BELOW . . .” printing is given and execution STOPs. In this case, "AT POSI- 
TION xxx” has no relevance. 

Subroutine date : September 12, 1969. 
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APPENDIX D 

LANGLEY LIBRARY SUBROUTINE DISCOT 


Language : FORTRAN 

Purpose : DISCOT performs single or double interpolation for continuous or discontinuous functions. 
Given a table of some function y with two independent variables, x and z, this subroutine performs 
K x th- and K z th-order interpolation to calculate the dependent variable. In this subroutine all single- 
line functions are read in as two separate arrays and all multiline functions are read in as three 
separate arrays; that is, 


x i 


L) 

y j 

(i = 1, 2, . . 

.,M) 

z k 

(k-1,2,. . 

.,N> 


Use : CALL DISCOT (XA, ZA, TABX, TABY, TABZ, NC, NY, NZ, ANS) 

XA The x argument 

ZA The z argument (may be the same name as x on single lines) 

TABX A one-dimensional array of x values 

TABY A one -dimensional array of y values 

TABZ A one-dimensional array of z values 

NC A control word that consists of a sign (+ or -) and three digits. The control word is formed 

as follows: 

(1) If NX = NY, the sign is negative. If NX * NY, then NX is computed by DISCOT as 
NX = NY/Nz, and the sign is positive and may be omitted if desired. 

(2) A one in the hundreds position of the word indicates that no extrapolation occurs above 
z max* With a zero in this position, extrapolation occurs when z > z max . The zero 
may be omitted if desired. 

(3) A digit (1 to 7) in the tens position of the word indicates the order of interpolation in 
the x-direction. 

(4) A digit (1 to 7) in the units position of the word indicates the order of interpolation in 
the z-direction. 

NY The number of points in y array 

NZ The number of points in z array 

ANS The dependent variable y 
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The following programs will illustrate various ways to use DISCOT: 

CASE I: Given y = f(x) 

NY = 50 

NX (number of points in x array) = NY 
Extrapolation when z > z max 
Second-order interpolation in x-direction 
No interpolation in z-direction 
Control word = -020 
DIMENSION TABX (50), TABY (50) 

I FORMAT (8E 9.5) 

READ (5,1) TABX, TABY 
READ (5,1) XA 

CALL DISCOT (XA, XA, TABX, TABY, TABY, -020, 50, 0, ANS) 

CASE II: Given y = f(x,z) 

NY = 800 
NZ = 10 

NX = NY/NZ (computed by DISCOT) 

Extrapolation when z > z max 
Linear interpolation in x-direction 
Linear interpolation in z-direction 
Control word = 11 

DIMENSION TABX (800), TABY (800), TABZ (10) 

1 FORMAT (8E 9.5) 

READ (5,1) TABX, TABY, TABZ 
READ (5,1) XA, ZA 

CALL DISCOT (XA, ZA, TABX, TABY, TABZ, 11, 800, 10, ANS) 

CASE III: Given y = f(x,z) 

NY = 800 
NZ = 10 
NX = NY 

Extrapolation when z > z max 
Seventh-order interpolation in x-direction 
Third-order interpolation in z-direction 
Control word = -73 

DIMENSION TABX (800), TABY (800), TABZ (10) 

1 FORMAT (8E 9.5) 

READ (5,1) TABX, TABY, TABZ 
READ (5,1) XA, ZA 

CALL DISCOT (XA, ZA, TABX, TABY, TABZ, -73, 800, 10, ANS) 

CASE IV: Same as Case in with no extrapolation above z max . Control word = -173 

CALL DISCOT (XA, ZA, TABX, TABY, TABZ, -173, 800, 10, ANS) 


41 



APPENDIX D - Continued 


Restrictions : See rule (5c) of section "Method" for restrictions on tabulating arrays and discontinuous 
functions. The order of interpolation in the x- and z-directions may be from 1 to 7. The following 
subprograms are used by DISCOT: UNS, DISSER, LAGRAN. 


Method : Lagrange’s interpolation formula is used in both the x- and z-directions for interpolation. This 
method is explained in detail in reference (a) of this subroutine. For a search in either the x- or 
z-direction, the following rules are observed: 

(1) If x < Xj, the routine chooses the following points for extrapolation: 

x i ,x 2 x k + i and yvh’ • ■ •» yk+i 

(2) If x > x n , the routine chooses the following points for extrapolation: 

x n-k’ x n-k+l x n y n -k>Vk + l y n 

(3) If x % x n , the routine chooses the following points for interpolation: 

When k is odd, 


x _k±!’ x 
2 



k + i and y 


i.sit+k 
2 


k+l ,y i k+1 |1 ’ 


When k is even, 






and 


y 

i- 





(4) If any of the subscripts in rule (3) become negative or greater than n (number of 
points), rules (1) and (2) apply. When discontinuous functions are tabulated, the indepen- 
dent variable at the point of discontinuity is repeated. 

(5) The subroutine will automatically examine the points selected before interpolation and if 
there is a discontinuity, the following rules apply. Let x d and x^ +1 be the point of 
discontinuity. 

(a) If x ^ x^, points previously chosen are modified for interpolation as shown: 

x d-k-*d-k+i' • • -’ x d yd-k>yd-k + i> • • **yd 

(b) If x > x d , points previously chosen are modified for interpolation as shown: 

x d+l* x d+2 x d+k and yd+l> y d+2> * * •» y d+k 

(c) When tabulating discontinuous functions, there must always be k + 1 points above 
and below the discontinuity in order to get proper interpolation. 

(6) When tabulating arrays for this subroutine, both independent variables must be in 
ascending order. 
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APPENDIX D - Concluded 


(7) In some engineering programs with many tables, it is quite desirable to read in one 
array of x values that could be used for all lines of a multiline function or different 
functions. Even though this situation is not always applicable, the subroutine has been 
written to handle it. This procedure not only saves much time in preparing tabular 
data, but also can save many locations previously used when every y coordinate had 
to have a corresponding x coordinate. Another additional feature that may be useful 
is the possibility of a multiline function with no extrapolation above the top line. 

Accuracy : A function of the order of interpolation used. 

Reference : (a) Nielsen, Kaj L.: Methods in Numerical Analysis. The Macmillan Co., c.1956. 

Storage : 555g locations. 

Subprograms used : UNS 40g locations. 

DISSER 110g locations. 

LAGRAN 55g locations. 

Subroutine date : August 1, 1968. 
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